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We investigate the dynamics of travelling oscillating solitons of the cubic NLS equa- 
tion under an external spatiotemporal forcing of the form f(x,t) = aexp[iK(t)x]. 
For the case of time-independent forcing a stability criterion for these solitons, which 
is based on a collective coordinate theory, was recently conjectured. We show that the 
proposed criterion has a limited applicability and present a refined criterion which is 
generally applicable, as confirmed by direct simulations. This includes more general 
situations where K (i) is harmonic or biharmonic, with or without a damping term in 
the NLS equation. The refined criterion states that the soliton will be unstable if the 
"stability curve" p(v), where p(t) and v(t) are the normalized momentum and the 
velocity of the soliton, has a section with a negative slope. Moreover, for the case of 
constant K and zero damping we use the collective coordinate solutions to compute a 
"phase portrait" of the soliton where its dynamics is represented by two-dimensional 
projections of its trajectories in the four-dimensional space of collective coordinates. 
We conjecture, and confirm by simulations, that the soliton is unstable if a section 
of the resulting closed curve on the portrait has a negative sense of rotation. 
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I. INTRODUCTION 

The externally driven, nonlinear Schrodinger (NLS) equation arises in many applications, 
for example charge density waves [1], long Josephson junctions [2], optical fibers [3-5] or 
plasmas driven by rf fields [6]. We use the NLS in the form 

iu t + u xx + 2|w| 2 -u + 5u = R[u(x, t); x, t], (1) 

with the perturbation 

R = f(x,t)-i/3u(x,t), (2) 

where f(x, t) is a direct (external) driving force and the term with j3 accounts for dissipation. 
Different forms of the driving force were considered: e.g., ac driving / = eexp(iwt) [1, 7, 8] 
or driving by a plane wave / = eexp[i(kx — out)] [5, 9]. Moreover, / = eexp[ig(x,t) — iut], 
where g is a function of x — vt, was considered, but no localized solutions were discussed [9]. 

The present paper continues the analysis [10] of the soliton dynamics under the spa- 
tiotemporal driving 

f{x,t) = ae iK M*. (3) 

A discrete version of Eq. (1) was used to model nonlinear optical waveguide arrays, in which 
discrete cavity solitons can be excited [11]. In that application, 5 is the cavity detuning 
parameter, and f(x, t) is replaced with f n (t) = a exp(i(p in n) , where n numbers the resonators 
and 4>in{t) is the incident angle of the laser pump light. A biharmonic function <pi n {t) was 
used in order to generate a ratchet effect [12] . In the present paper we also obtain a ratchet 
effect by using a biharmonic driving (Section V). This is interesting because there are only a 
few reports on ratchets with nontopological solitons [12-14]; most of the literature concerns 
ratchets with topological solitons, e.g [15-19]. 

In Ref. [10], Eq. (1) was simulated using the 1-soliton solution of the unperturbed NLS 
as the initial condition. The soliton's position, velocity, amplitude and phase served as 
parameters of the initial condition (IC). 

In the case of zero damping and time- independent, spatially periodic driving of the form 
f(x) = exp(iKx) the resulting solitons were observed [10] to display periodic oscillations 
of their positions, velocities, amplitudes and phases. Although the driving force has zero 
spatial average, the soliton's net motion is unidirectional. (This contrasts with the case of 
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a perturbation V(x)u with periodic V(x), where the soliton performs oscillations about a 
minimum of V(x) [20]). 

A large number of sample points in the parameter space (a, K, 5) were examined by 
varying the initial amplitude i]q, with the other initial conditions kept fixed. The initial 
configuration was seen to evolve into a stable soliton only when t]q was taken from one of 
the "stability windows" . For i] outside the stability windows, the solitonic initial condition 
was observed to decay or break into two or more fragments which would subsequently decay 
[10]. 

As a first step towards understanding the observed dynamics of solitons, the authors 
of Ref. [10] proposed an empirical stability criterion based on a collective coordinate (CC) 
description. The collective coordinates analysis produces a set of coupled nonlinear ODEs for 
the soliton's position q, amplitude n, normalized momentum p and phase $. An approximate 
solution of this dynamical system is given by trigonometric functions and can be obtained 
explicitly, except when the initial condition i] is chosen near one of the stability boundaries. 
In the latter case the collective coordinates equations had to be analysed numerically and 
their solutions were found to be highly anharmonic. 

We have positively tested the predictions of the proposed stability criterion by simulations 
(numerical solutions of the full NLS Eq. (1)) for many classes of initial conditions. However, 
the tests are negative, if the initial momentum is too large, i.e. p > K for positive K (or 
Po < K for negative K). In this paper we therefore conjecture a refined stability criterion, 
which we show makes correct predictions not only for the case K = constant (with all 
classes of initial conditions), but also for harmonic and biharmonic K(t). The new criterion 
is a sufficient condition which states that the soliton will be unstable in simulations, if 
the "stability curve" p(y) has a branch with negative slope. This curve is obtained as a 
parametric plot of the normalized momentum 

P=^ (4) 

p N(ty { ) 

of the soliton versus its velocity, 

v = ?(*)■ (5) 
Here P = 4np is the canonical momentum of the soliton, N = An is the norm which is 
canonically conjugated to the soliton's phase <&(£), see Section II. In the old criterion [10] 
the stability curve was defined as P(v). At the end of Section III we present an example 
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which demonstrates by an analytical calculation that the normalized momentum p, instead 
of the canonical momentum P, has to be used for the stability criterion. 

It is important to emphasize that the soliton's stability or instability is judged not on 
the basis of the stability of solutions to the collective coordinates equations. (The latter are 
stable in most cases). The soliton's stability is rather decided on the basis of some of its 
properties which are captured by the p(v) curve of the corresponding collective coordinates 
solutions. The proposed empirical criterion reproduces the numerically observed positions 
of the stability windows to an accuracy of better than 1%, despite the complexity of the 
stability diagram in the parameter space [10]. 

The stability criteria for the homogeneous (translation invariant) NLS equation available 
in the literature are restricted to (a) bright solitons, i.e. solutions decaying to zero at the 
spatial infinities, with time dependencies of the form e iAi (and those reducible to this form by 
a Galileian transformation); (b) traveling dark solitons, i.e. solutions approaching nonzero 
constant values as x — > ±oo. The criteria are insensitive to the particular form of the 
nonlinearity as long as it is conservative and [/(l)-invariant, i.e. as long as the equation 
does not include any damping or driving terms. 

In the case of the bright solitons of the form u(x,t) = u s (x)e tAt , the Vakhitov-Kolokolov 
criterion states that if the corresponding energy Hessian has only one negative eigenvalue, 
then the soliton is stable if dN/dA > and unstable otherwise [21-23]. Here N = J \u\ 2 dx; 
depending on the physical context, N is referred to as the number of particles contained in 
the soliton or the total power of the optical beam. (See also [24] for the energy- versus-number 
of particles formulation of this criterion.) 

In the case of dark solitons of the form u(x,t) = u(x — Vt), with \u\ 2 — > p as \x\ — > oo, 
a similar criterion [25-28] involves the (renormalised) field momentum, 

p = \J («;« - «*«•) (i - dx. (6) 

The dark soliton travelling at the constant velocity V is stable if dP/dV < and unstable 
otherwise. 

Some parts of the stability analysis of the travelling dark solitons [27] can be carried over 
to the case of the travelling solitons of the NLS with a driving term. Namely, one can show 
[29] that a linearised eigenvalue crosses from the negative to the positive real axis at the 
value V where dP/dV = 0. The sign of the derivative dP/dV required for stability depends 
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on the type of the soliton; some classes of solitons require dP/dV < 0, whereas other classes 
are stable when dP/dV > 0. (An additional complication is the presence of oscillatory 
instabilities where two eigenvalues collide on the imaginary axis and acquire opposite real 
parts. The oscillatory instabilities do not affect the sign of dP/dV.) 

In these analyses, each point of the curve P(V) represents a soliton traveling at a par- 
ticular constant velocity V; therefore the curve is a characteristic of the whole family of 
solitons. The values of V where dP/dV = break the family into parts with different sta- 
bility properties. In contrast to this, each oscillatory solution of the collective coordinates 
equations [10] has its own, individual, p(v)-curve the whole of which is traced periodically 
in time. The shape of this curve determines whether the corresponding soliton is stable or 
not. 

The present paper has several goals: First, we propose a refined stability criterion. Sec- 
ond, we study the internal structure of the instability regions. We will demonstrate that 
these regions consist of subregions characterized by instabilities of different types. The ex- 
istence of the subregions will be predicted by the analysis of the reduced dynamical system 
and confirmed by direct simulations of the full PDE (Section III). In obtaining the reduced 
dynamical system we modify the original collective coordinates approach of Ref. [10] (Sec- 
tion II). In addition to producing bounded trajectories (a property essential for the stability 
analysis), the modified approach provides a much easier derivation of the canonical soliton 
momentum and the Hamilton function in terms of the canonical variables (Section II). 

Third, we demonstrate that a certain "phase portrait" of the soliton on the complex 
plane can be used as an alternative stability diagnostic (Section III). However, the phase 
portrait requires the phase of the soliton to be periodic in time. This can be achieved by the 
above mentioned modification of the original collective coordinates approach [10] in which 
the phase was not periodic, in contrast to the other three collective coordinates. 

Finally, we explore the applicability of our refined stability criterion to inhomogeneous 
forcings of the form f(x,t) = aexp(iK(t)x) in Eq. (3). We will start with a harmonically 
varying K(t), with and without the damping term in the right-hand side of (1) (Section 
IV). After that, in Section V, we will consider a biharmonic K(t) with a broken temporal 
symmetry. (The temporal symmetry breaking will accompany the breaking of the spatial 
symmetry by the inhomogeneous driving.) 
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II. MODIFIED COLLECTIVE COORDINATE THEORY 

The one-soliton solution of the unperturbed NLS is given by [30] 

u(x, t) = 2i V sech[2r/(:r - C)]e" i(2fx+0) , (7) 

where rj and £ are real parameters (r/ > 0); C = Co — 4£t gives the coordinate of the soliton's 
center, and = 0o + (4£ 2 — 4r/ 2 — 5)t is the soliton's phase. The collective coordinates theory 
of Ref. [10] assumed that for sufficiently small perturbations R in Eq. (2) the soliton shape 
and dynamics can be described, approximately, by Eq. (7), where r/(t), ((t) and (f)(t) 
are functions of time. 

We now show that the following modification of this ansatz [31, 32] provides a considerable 
improvement of the collective coordinates theory of Ref. [10]: 

u(x, t) = 2i V sech[2r](x - q)]e l[p(x - q) -' s ' ] , (8) 

by setting — 2£ = p, ( = q, and = $ — 2£( = $ + pq. Here only the last replacement 
is essential for the above mentioned improvement of the collective coordinates theory. The 
four collective coordinates equations of Ref. [10] are replaced with 

71 

f] = —2(3i] — —a sechA cos B, (9) 

7T 2 O, 

q = 2p-\ -sechA tanh A sin B, (10) 

8 r] 2 

p = -2aA sechA cos B, (11) 

$ + pq = p 2 - Arf - 5 H sechA tanh A sin B, (12) 

2 rj 



with 



A(t) = ^ f) [K(t)-p(t)l (13) 
B(t)=$(t) + K(t)q(t). (14) 



The new formulation has the following advantages 
1. Consider the Lagrangian for Eqs. (9)-(12): 
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L = 4r/$ + 4r]pq - 4r]p 2 + — r/ 3 + A5r) - 2ixa sechA sin B. (15) 



The momentum conjugate to the phase $ is 

dL 



,, H (16) 
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which is equal to the norm ( f \u\ 2 dx) of the waveform (8). The momentum conjugate to 
the soliton's position is 

dL , x 

— = 4 V p. (17) 

The advantage of the new formulation is that this is equal to the field momentum of the 
configuration (8) 



P — - J {u* x u — u x u*)dx, (18) 
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whereas in Ref. [10] the second canonical momentum was defined by dL/dp = —Arjq which 
did not have any obvious physical interpretation. 

If the dissipative term —if3u in (2) has a nonzero coefficient, we have to use the generalised 
Euler-Lagrange formalism with the dissipation function 

/+oo 
(uu* t - u*u t )dx. (19) 
-oo 

Substituting (8) in (19), we obtain 

F = -8(3ri(<5>+pq). (20) 

The generalised Euler-Lagrange equations are 

d_dL_ _ dl^ _ dF_ 
dt dif ~ d4> ' 

where \& represents each of the four collective coordinates 77, q, p, and <E>. 

2. Since P = 4rjp is the canonically conjugate momentum for q, the Legendre transform 
to the canonical Hamiltonian is easily performed: H = N<& + Pq — L. This gives 

H = ^-P 2 - — N 3 -5N + 2nasechA sin B. (22) 
12 v ' 

In Ref. [10] this Hamiltonian could only be obtained via a canonical transformation. 

When the forcing f(x) = ae tKx is time- independent and damping /3 — 0, the collective 
coordinates rj and p perform periodic oscillations, whereas q(t) and are given by periodic 
functions superimposed over linearly growing functions of t. In contrast to this, the variable 
= $ + pq used in Ref. [10] will obviously exhibit oscillations with a (linearly) growing 
amplitude. 
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III. TIME-INDEPENDENT, SPATIALLY PERIODIC FORCE 

Our approach consists in the numerical solution of the collective coordinates equations 
(9)-(12) for representative values of the parameters and the initial conditions i]q, q , p and 
$o- Each collective coordinates orbit is then used to compute p(t) and v(i) in Eqs. (4)- (5) 
and to plot p against v. If some part of this "stability curve" has a negative slope, we predict 
that the soliton will become unstable in simulations of the PDE (1) starting with the initial 
condition (8), with the same rj , p , q and $ . 

Since the collective coordinates approximation can work only for small perturbations, we 
choose a small driving amplitude a = 0.05 for f(x) = ae lKx . We also take K = 0.1, which 
means that the spatial period of the forcing L — 2n/K >> 1. We are interested in periodic 
solutions and therefore we set the damping parameter (3 = 0. Damped oscillatory solutions 
were already considered in Ref. [10]. 

For 5 > 0, the p(v) curve predicts only unstable solitons which is confirmed by the 
simulations. For 5 < 0, there are typically several stability regions which grow as \S\ is 
increased, while the parameters a and K are fixed [10]. We concentrate here on the simplest 
case with only one stability and one instability region. Namely, we choose 5 = — 1, q = 
p = $o = for which the soliton solutions are predicted to be stable if 7] > r]P = 0.684 
(Fig. la) and unstable for rjo < rj^K This is confirmed by our simulations of the PDE (1) to 
an accuracy of better than 1% in rj^K 

The range of initial amplitudes i]q for which solutions of Eqs. (9)-(12) feature a p(v) 

(2) 

curve with a descending branch, can be divided into two subintervals, < 770 < ffc and 
f]c" > < rj < r]i\ where r/^ = 0.288. In the upper subinterval, (r)]?\r]^), the curve p(v) 
exhibits two long branches, one with a positive and the other with a negative slope (Fig. lb). 
In the lower subinterval, the positive-slope branch is short and very steep (Fig. lc). The 
difference in the shape of the stability curves suggests different types of instability in the 
two subregions; however, in which way the instabilities are different cannot be deduced from 
the p(v) curve alone. To gain further insight into this difference, we plot a phase portrait for 
the dynamical system (9)-(12). The vertical and horizontal axes in the portrait are chosen 
so that they admit a clear interpretation in terms of the full PDE, Eq. (1). To this end, we 
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FIG. 1. Stability curve p(y) corresponding to a = 0.05, K = 0.1, S = — 1 and /3 = 0. The initial 
conditions for the Eqs. (9)-(12) were qo = po = $o = and a) 770 = 0.8, b) 770 = 0.65, c) 770 = 0.1. 
The integration time tf = 1000. 

first transform to the frame of reference moving with the velocity Vf. 

u(x, t) = t)e iKx , X = x - V f t, (23) 

where Vf = 2K. Eq. (1) is taken to be an NLS driven by a space-time independent external 
force: 

iVt + ^xx + 2|^| 2 ^ + = a. (24) 
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(This equation was previously studied in a different context [7, 8, 33, 34] and two static 
soliton solutions were obtained explicitly [34]. Unlike [7, 8, 33, 34], we focus here on moving 
solitons of Eq. (24).) 

Under the transformation (23), the collective coordinates ansatz (8) becomes 

*(X,t) = 2ir)sech[2r)(X + V f t - 5 )] e -*[(*-p)(^/0+w+*]. ( 2 5) 

The rj(t) and p(i) components of the oscillatory solutions of (9)-(12) are periodic with period 
T, whereas q and $ are of the form q(t) = vt + q (p) (t), ®(t) = -at + $ (p) (t), where q (p \t) 
and <E>( p )(t) are T-periodic functions and a is a constant [10]. The corresponding soliton (8), 

(25) has the mean velocity v in the original frame of reference, and v — Vf in the moving 
frame. 

At the point x = vt [or, equivalently, at X = (v — Vf)t], the function (25) has the following 
time dependence: 

= 2irj sech[2r](vt - g )] e -*[***-p(™-«)+*]. ^ 

The function (26) is a collective coordinates counterpart of the ^ field at the centre of the 
soliton solution of Eq. (24). Comparing Eq. (26) to the function \1/(X, t)\x=n-v f t obtained in 
the direct numerical simulations of the full PDE (1), one can assess the validity and accuracy 
of the collective coordinates approximation. For this reason, we choose the complex function 

(26) as a representative of the four-dimensional dynamics, and plot its real versus imaginary 
part to generate the corresponding phase portrait. The soliton dynamics is described by the 
resulting orbits of the phase portrait. 

One can readily verify that these orbits are closed. Indeed, the modulus of the function 
(26) is periodic with period T. Therefore, to demonstrate the closure, one just needs to 
show that the argument of \1/ changes by an integer multiple of 2tc over the period. We 
have arg\I> = (a — Kv)t — p^q^ — where the last two terms are T-periodic. As for 
the first term, the constants v and a are found as the coefficients of the linearly-growing 
components of q(t) and $(£), respectively. The numerical solution of Eqs. (9)-(12) verifies 
that a — Kv = 2n/T in the stability range r]o > rftt\ that a — Kv = —2n/T in the lower 
instability subinterval < r] < rj c , and Kv — a = in the upper instability subinterval 
Vc 2 ^ < Vo < Vc 1 ^ ■ These relations between v and a hold to a numerical accuracy of O(10~ 5 ). 

Trajectories resulting from initial conditions in the interval i] > i]^ are ellipses, with a 
positive sense of rotation (Fig. 2). The ellipses enclose a stable and an unstable fixed point 
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on the real axis at about +1 and —1, respectively. (For the definition and calculation of 
these points see Appendix A). Fig. 3a compares the soliton amplitude rj(t) from collective 
coordinates theory to the amplitude measured in the direct simulations of Eq. (1). 




FIG. 2. The phase portrait of the system (9)-(12) with a, K, S and (3 as in Fig. (1). Shown is 
lm^(X = vt — Vft,t) versus Re^(X = vt — Vft,t). The large ellipse corresponds to rjo = 0.8, the 
horseshoe to t]q = 0.65 and the small ellipse to t]q = 0.1. Other initial conditions are as in Fig. (1). 
The separatrix is shown by the dotted curve. The filled and open circles are stable and unstable 
fixed points, respectively. 

For the upper instability subinterval 7]^ < r] < rjc the phase trajectory is a horseshoe 
(Fig. 2). This curve consists of an outer part with a positive sense of rotation and an inner 
part with a negative sense of rotation relative to the origin. The two parts are correlated 
with the two branches with positive and negative slopes, respectively, of the p(v)-curve in 
Fig. lb. The soliton instability is seen in the simulation result in Fig. 3b. Note that the first 
harmonic vanishes after about 30 time units, while the second harmonic persists. Eventually 
the soliton decays: the amplitude approaches zero while the width tends to infinity. 



For the lower instability interval < rj < rjc the situation is quite different, both in the 
collective coordinates theory and in the simulations: The phase portrait features an ellipse, 
but with the negative sense of rotation (Fig. 2). Moreover, the ellipse is much smaller than 
the one arising in the stability region so that it encloses only one fixed point. This fits with 
the simulations in which the soliton remains metastable for a relatively long time, exhibiting 
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FIG. 3. Soliton amplitude r)(t) from the collective coordinates theory (solid lines) and from the 
simulations (dashed lines). The parameters and the initial conditions are the same as in Fig. 1. a) 
7/o = 0.8, b) 7/o = 0.65, c) 7/o = 0.1 (shown are results for early times < t < 100); d) 7/0 = 0.1 
(shown are simulation results for late times 600 < t < 1000.) 
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FIG. 4. Stability diagram near r/c = 0.684, with q = <J? = 0. Parameters: a = 0.05, K = 0.1, 
5 = —1, /3 = 0. Circles: unstable soliton. Plus: stable soliton. Dashed line: for po = ^\ the p(v) 
curve is a point. 

a periodic modulation of the oscillation amplitude (Fig. 3c), but then the instability sets in 
(Fig. 3d). 

So far we have varied 770, with p = q = $ = kept fixed. We now consider the stability 
diagram in the 770-po plane near the critical value rfc \ which separates the stability interval 
from the upper instability interval (above). Fig. 4 shows that a finite value of the normalized 
momentum p stabilizes the soliton and therefore the stability region is enlarged. The curve 
which separates stability and instability regions is roughly a parabola. 

Finally we would like to emphasize the crucial importance of using the normalized mo- 
mentum p in the stability analysis — rather than the canonical momentum P (as was pro- 
posed in [10]). We have established that the empirical stability criterion suggested in [10] 
disagrees with the results of numerical simulations when the initial normalized momentum 
is too large, i.e. p > K (for positive K). Let, for instance, the parameters of the equation 
take the same values as in Fig. 1 (a = 0.05, K = 0.1, 5 = — 1 and (5 = 0), and take the same 
initial conditions as for the stable stationary solution in the Appendix (770 = 0.5a/K 2 — 5, 
$0 = tt/2, qo = 0), except that this time po = K + d, where < d < 0.2. In this case 
the numerical solutions of the collective coordinates equations can be represented in a very 
good approximation by p(t) = po + a p (l — cos fit) and v(t) = v + a v (l — cos fit), where po, 
v , a p , a v > 0. Thus p(v) is a straight line with slope a p /a v > 0. This predicts stability, the 
same as the orbit in the phase portrait which is a small ellipse with positive sense of rotation 
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around the stable fixed point at about +1 on the real axis. The stability is confirmed by 
simulations. However, when the momentum P = 4r)p is used the situation is different. r)(t) 
can be expressed viap(t) by using the exact relation r\ = i] (p — K) / (p — K) {p ^ K), which 
is obtained from Eqs. (9) and (11) where one integration has been carried out. Finally, p(v) 
from above is inserted and one can see that P decreases when v increases and viceversa. 
Thus the slope dP/dv < predicts instability which disagrees with the simulations. 

IV. HARMONIC K(t) 

As stated in the Introduction, one of the aims of this paper is to verify whether the sta- 
bility criterion p'(v) > remains applicable to time- dependent forces of the form f(x,t) = 
ae iK{t)x _ j n j-jj-g sec ti on we consider the case of a harmonically modulated forcing wavenum- 
ber: 

K(t) = ksm(cot + 6), (27) 

first without a damping term in the NLS ((3 = 0), then with the damping (/? > 0). 

We choose the same parameters as in Section III: a = 0.05, k = —0.1 which implies 
| if | <C 1. In order to be in the adiabatic regime we choose a small modulation frequency 
u = 0.02. Finally, we let 6 = and choose initial conditions q = Po = 0, rj = 1, and 
$ = tt/2. 

The numerical solutions of the collective coordinates equations (9)-(12) exhibit oscilla- 
tions with three very different frequencies in their spectrum. This is most explicit in the 
behaviour of q(t) (Fig. 5 a,b). First, there are intrinsic oscillations with the frequency uf, 
these have a period Tj of the order of 10, similarly to the oscillations in the case of the 
constant K discussed in the previous section. Second, there are oscillations with the driving 
frequency u whose period = 2tt/u 314. Finally, there are oscillations with a very low 
frequency ui and very long period TJ 8000 (Fig. 5b). The resulting stability curve p(v) 
exhibits many small loops which have a short section with a negative slope. An example 
is given in Fig. 5c. (For clarity the curve is plotted only over a short time interval). The 
negative slope predicts instability; this is confirmed by our simulations of the full PDE. 

Stable solitons can be obtained by changing r/ in such a way that the loops do not arise. 
This is achieved by suppressing the intrinsic oscillations, since their period T { w 10 is of the 
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FIG. 5. Collective coordinates results for harmonic K(t), no damping, a = 0.05, k = —0.1, 
co = 0.02, = 0, S = —3, j3 = 0, qo = po = 0, $o = tt/2, % = 1. a) q(t) exhibits Wj-oscillations 
modulated by the frequency oj. Shown is the interval < t < Td = 2ir/u. b) q(t) exhibits oj- 
oscillations modulated by the frequency uji = 2n jT\. Here < t < 8000. c) Stability curve p(v) for 
— 5 < t < Tj/2 + 10. The arrow points to the section of the curve with a negative slope loop. 
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same order as the time scale of the loops, see Fig. 5c. The intrinsic oscillations disappear 
when we choose rj = (Fig. 6a,b). In this case r](t) performs very small oscillations 

around r] and the two dominant terms on the r.h.s. of Eq. (12), namely — 4t]q and —8, 
cancel each other. Fig. 6c demonstrates that the small loops have indeed disappeared. 

The resolution of Fig. 6c does not allow verification of whether there are sections with 
negative slope near the turning points of the stability curve. We now show that there cannot 
be any, as the curve develops cusps at the turning points. Consider the region around one 
of the maxima (or minima) of the w-oscillations of the collective coordinates (Fig. 6a). The 
functions q(t), r)(t) etc. are not symmetric with respect to t m (position of the extremum) 
due to the existence of the very slow o^-oscillations. The same holds for p; hence 



with Ci ^ C r . For the velocity v(t) = q(t) the asymmetry is negligible, compared to the 
asymmetry of pit), because the time derivative q contains a factor uji <C 1. Thus v(t) = 
v m — b{t — t m ) 2 for both t <t m and t > t m . Eliminating t we obtain p = p m — Q(v m — v)/b 
as v increases up to its maximal value v m , and v = v m — C r {y m — v)/b as v decreases from 
v m . Thus the stability curve p{y) has a cusp with two different (but positive) slopes C\jb 
and C r /b at the turning point v = v m . The absence of segments of the curve with dp/dv < 
predicts stability for the soliton. This is confirmed by the simulations of the PDE (1). 

When the damping term — ifiu is included in the r.h.s. of the NLS equation (1), the 
collective coordinates dynamics simplifies. Namely, both the intrinsic oscillations and the 
low-frequency oscillations are damped out from solutions of the collective coordinates equa- 
tions after a transient time t tr = 1/(5. After this transient, all collective coordinates oscilla- 
tions become locked to the driving frequency u. The stability curve in this case consists of 
two nearly-straight lines which form sharp cusps at both ends (Fig. 7a). Thus there are no 
sections with a negative slope and the soliton is predicted to be stable. This is confirmed by 
the simulations of the PDE. For long times (t ^> t tr ) the average soliton velocity v slowly 
approaches zero (Fig. 7b); this behaviour is independent of the initial conditions. Thus there 
is no unidirectional motion of the soliton for long times; the reason will be established in 
the next section. 




Pm - Q(t - t m f for t < t 
p m - C r (t - t m ) 2 for t > t. 



in 



m 



(28) 
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0.8663 




time 




time 




V 



FIG. 6. Collective coordinates results for harmonic K(t), no damping. Same parameters and initial 
conditions as in Fig. 5, but 770 = y/—S/2. (For this choice the intrinsic Wj-oscillations vanish.) a) 
The amplitude r](t) exhibits no ^-oscillations, only w-oscillations modulated by ^/-oscillations 
which are hardly visible for < t < 1000. b) For rj(t) for < t < 8000 both oj- and w^-oscillations 
are visible, c) Stability curve for < t < T^. 



As (3 is decreased, the stability curve becomes wider and the decay of v to zero faster. 
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On the contrary, as j3 is increased, the stability curve becomes narrower, and v decreases to 
zero more slowly. However, for (5 above a critical value f3 c (f3 c « 0.035 for the parameters of 
Fig. 7) the collective coordinates solutions become unstable. Direct simulations also confirm 
the soliton's instability. 





FIG. 7. Collective coordinates results for harmonic K(t), nonzero damping. Same data as in Fig. 
6, but T] = 0.866 « y/=8/2, /3 = 0.01. a) stability curve p(v) for 3T d < t < 4T d . b) q(t) for 
10000 < t < 100000. 



V. BIHARMONIC DRIVING: RATCHETS 



The simplest ratchet models consider a pointdike particle in a periodic potential driven 
by an AC force f(t). Under certain conditions related to the breaking of symmetries, unidi- 
rectional motion of the particle can take place despite the applied force having zero temporal 
average [35-39]. Particle ratchets were generalized to nonlinear field theoretic systems, in 
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which particles are replaced by solitons [40-46]. In particular, solitons in nonlinear Klein- 
Gordon systems can move on the average in one direction, if either a temporal or a spatial 
symmetry is broken. 

A temporal symmetry, namely a time-shift symmetry, is broken by a biharmonic force [47, 
48]. In this case the mechanism of the ratchet effect was clarified by a collective coordinates 
theory employing the soliton position and width as collective coordinates [15, 16, 49]. Due 
to the coupling between the translational and internal degrees of freedom, energy is pumped 
nonuniformly into the system, generating a unidirectional motion. The breaking of the 
time-shift symmetry gives rise to a resonance mechanism that is present whenever the soliton 
oscillation spectrum comprises at least one of the frequency components of the driving force. 

In this section we investigate whether the NLS solitons show a behavior similar to 
the Klein-Gordon kinks. This is particularly interesting since the NLS solitons are non- 
topological, whereas the vast majority of reports on soliton ratchets have so far focussed on 
topological solitons. 

We consider the NLS (1) where the perturbation 

R = f(x,t) -ipu, (29) 



has the form 



with 



f(x, t) = ai e iKlit)x + a 2 e iK2{t)x , (30) 



K l = fcx sin(wt), K 2 = k 2 sm(2cot + 9). (31) 

Consider first the single-harmonic case: a 2 = 0. When t 3> t tr = 1//3, the soliton oscillations 
are locked to the driving frequency uj and are independent of the initial conditions (see 
Section IV). Thus there exists a global solitonic attractor. 

We now perform a symmetry analysis [12, 48]. The perturbed NLS is invariant under the 
symmetry operation 

S : t^t + T/2, x i — y -x. (32) 

At the same time, the transformation S changes the sign of the soliton velocity v(t) = X(t). 
The soliton position is defined by 

f +0 ° dx xp(x, t) 

x{t) = 7 (33) 

J-oo dxp{x,t) 
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with 

p(x,t) = \\u(x,t)\ 2 -\u bg (x,t)\ 2 \. (34) 

Here Ub g (x,t) = a^it) exp^iKxx) is the background field to which the soliton decays as 
\x\ — > oo [10, 12]. When \u(x, t)\ 2 from the simulations is plotted, the soliton sits on a shelf 
with homogeneous intensity \ab g (t)\ 2 . The shelf has little influence on the soliton dynamics 
[10]; this is why the collective coordinates theory is in a good agreement with simulations, 
despite ignoring the presence of the background. 

Since the attractor is global, the transformation S maps it onto itself. This implies that 
the average velocity on the attractor is zero. The soliton performs periodic oscillations about 
its equilibrium position which are reproduced by the collective coordinates theory (Fig. 7b). 

In order to construct a solitonic ratchet, i.e. obtain a stable soliton with v ^ 0, it is 
necessary to break the invariance under the operation S. The simplest way to do this is 
to employ the biharmonic driving in Eq. (30) with a± ^ and a 2 ^ 0. The collective 
coordinates equations (9)-(12) can easily be extended to the case of the forcing function / 
including two terms. In particular, Eq. (9) is replaced with 

2 

r) = — 2/??7 — di — sechAj cos Bj, (35) 

i=i 

where 

4 = 4^[Ki(*)-p(*)l, (36) 
B i = Q(t) + K i (t)q(t), (37) 

while the Ki are as in Eqs. (31). The collective coordinates equations for q, p and $ are 
modified in a similar way. 

Since the collective coordinates description is accurate only for small perturbations, we 
take small driving amplitudes ai = a 2 = 0.05. We choose a very small driving frequency 
oo = 0.002 in order to remain in the adiabatic regime. If the damping coefficient /3 is chosen 
too large, the soliton amplitude i] quickly relaxes to zero while q(t) and <&(£) rapidly go to 
infinity. For example, for the parameters 5 = —3, k± — k 2 — k — 0.001, 9 = and the IC 
r] — 1, q — p — 0, $o = 7r/2, this instability occurs when f3 > 0.065. On the other hand, if 
f3 is chosen too small (e.g., f3 = 0.01), the average soliton velocity grows without bound over 
sufficiently long integration times (tf ~ 10 5 ). Thus we can expect a stable ratchet effect 
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only for intermediate values of f3, for instance (3 = 0.04. As we find that v ~ k, a larger 
ratchet effect can be obtained by increasing k. However, when k exceeds a certain critical 
value k c , the average velocity starts to grow slowly with time. (For the chosen parameter 
values, k c = 0.002). 

Using the parameter values and initial conditions for which the collective coordinates 
equations exhibit stable solutions we perform direct simulations of Eq. (1). Our aim is to 
test whether an initial waveform (8) will evolve into a stable solitary wave over the time 
t tr = 1/(5. However, it turns out that the initial structure evolves rapidly immediately after 
the start of the simulation and quickly decays to zero. 

In order to obtain stable solitary waves we need to improve the initial conditions. This 
can be achieved by setting the initial conditions equal to the mean values about which the 
collective coordinates oscillate, once the transients have elapsed. These mean values can be 
obtained from an approximate analytical solution of the collective coordinates equations: 
We let 

q = vt + C q , (38) 
V = P + Q, 

7] = fj + C v , 

$ = $ + c$, 

where C x are oscillations with amplitude a x and zero mean. We choose u = O(10~ 3 ), and 
k\ = k 2 = k = O(10 -3 ). The other parameters (5, a-y, a 2 , /?) do not have to be small for the 
following perturbation analysis and can therefore be chosen in 0(1). Substituting in the 
collective coordinates equations we retain only the leading terms in the perturbation series. 
This gives 

P = 0, p= a i*i(*) + ^(t) =0(10-3), (39) 
(ai + a 2 ) 

n = \^, s = o(io- 6 ), (40) 

$ = arccos ( - ) , a* = O(10" 3 ), (41) 

V7r(ai + a 2 ) J 

k 2ai . . k 2a 2 . . 

q = vt+ — (1 -cos(wt)) + — (42) 

to a\ + a 2 2tu ai + a 2 

x(l-cos(2wt + 0)), v = O(10~ 6 ). (43) 
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Eqs. (39)- (42) are in very good agreement with the numerical solution of the collective 
coordinates equations. We note that the constants p, fj and $ do not depend on the relative 
phase 9, but the variable components of q(t) and p(t) do. 

The improved initial conditions are now: p = p = 0, r) = fj, $ — Qo — 0. After a 
transient time the numerical trajectory settles to the solution (39)-(42). This yields v(t) = 
q = 2p and thus p(v) = |v is a straight line with positive slope. Our stability criterion 
predicts the stability of the soliton and simulations of the PDE confirm this (Fig. 8a,b). As 
v is very small, the ratchet effect is not visible on the time scale of Fig. 8a,b, but can be 
observed over longer simulation times tf « 30T d (Fig. 8c). The soliton collective coordinates 
amplitude rj(t) oscillates about 0.866025 in the interval [0.866024,0.866026], in agreement 
with Eq. (40). In the simulations, rj(t) oscillates about 0.87363 in [0.87305,0.87390]. 

The ratchet effect is also observed for higher driving frequencies (e.g., cu = 0.01 and 0.02). 
However, the average velocity v decreases in proportion to 1/ou, similar to the last two terms 
in Eq. (42)). It is important to emphasize here that we only succeeded in determining 
initial conditions for a stable soliton thanks to the availability of the explicit solution of 
the collective coordinates equations and our stability criterion. It would be very difficult to 
identify the corresponding small basin of attraction via numerical simulations of the PDE. 

Finally, we discuss the dependence of the average velocity v on the relative phase 9 in the 
biharmonic driving force (31). As expected for a ratchet system with biharmonic driving 
[48, 49], v(9) is sinusoidal with the period 2n. It attains its maximum value near 9 = and 
its small negative minimum value near 9 = ir. The size and shape of the basin of attraction 
around (fj, <5) also depend strongly on 9; this effect will be examined in a future work. 

VI. SUMMARY 

We have formulated a refined empirical stability criterion for the driven NLS solitons. 
Unlike stability criteria available in the literature, the new criterion is based on a Collective 
Coordinate (CC) description. Solving (analytically or numerically) evolution equations for 
the four collectives coordinates, we use the resulting trajectories to evaluate the normalized 
soliton momentum p(t) and the soliton velocity v(t). These give a parametric " stability 
curve" , p(v). 
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2000 4000 6000 8000 



time 





FIG. 8. Soliton position q(t): (a) from the collective coordinates theory, (b) from simulations, (c) 
the ratchet effect in the collective coordinates theory visible over long times (tf « 30 periods). 
Parameters: oi = a 2 = a = 0.05, h = k 2 = k = 0.002, 5 = -3, u = 0.002, = 0,0 = 0.08, with 
initial conditions t]q = fj, qo = po = 0, <&o = 



Whenever the curve p(v) has a section with a negative slope (dp/dv < 0), we observe the 
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instability of the soliton in direct numerical simulations. We therefore conjecture that the 
availability of a section with a negative slope is a sufficient condition for the instability of 
the soliton. We do not have a mathematical proof of this conjectured criterion; however we 
have verified it in a variety of situations using constant, harmonic and biharmonic functions 
K(t), with or without the damping term. 

Establishing a theoretical justification of this conjecture is a subject of future work, first 
for the cases of K(t) which we have considered in this paper, then for a general function 
K(t). One of the foreseen difficulties is related to the fact that the soliton solution of the 
driven NLS does not vanish as x — > ±00, because the perturbation f(x,t) = aexp[iK(t)x] 
does not decay to zero in these limits. On the other hand, our collective coordinates theory 
is based on a soliton ansatz which vanishes as x — > ±00. 

For the case of constant K and zero damping, all collective coordinates perform periodic 
motions. This allowed us to compute a phase portrait which consists of closed orbits on a 
complex plane. The soliton evolution is described by motion along one of these orbits. We 
observe that the sense of rotation of the orbit is correlated with the stability/instability of 
the soliton (determined in simulations of the full PDE). Namely, 

1. If the orbit is an ellipse with a positive sense of rotation, the soliton is stable. 

2. If the orbit is a horseshoe where the inner part has negative and the outer part positive 
sense of rotation, the soliton is unstable and desintegrates very quickly. 

3. If the orbit is an ellipse with a negative sense of rotation, the soliton remains metastable 
for a relatively long time but eventually desintegrates. 

An interesting question is whether our collective coordinate approach and our stability 
criteria can also be applied to NLS equations with a more general form of the nonlinearity. 
Work is in progress regarding the case of a nonlinearity with arbitrary power, {u*u) K , where 
k — 1 corresponds to the NLS of this paper. The unperturbed NLS has stable solitons for 
< k < 2 and it will be interesting to determine how the stability of the solitons is affected 
by the perturbation f(x,t). 
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VIII. APPENDIX A: FIXED POINTS OF THE PHASE PORTRAIT 

For the case of the time-independent force f(x) = ae tKx and zero damping, we adopt the 
following ansatz for stationary solutions of the collective coordinates equations: 

q(t) = v s t, 77(f) = r] s , pit) = p s , $(t) = $ s - a s t. (44) 
Eq. (9) yields cos B = which results in 

Kv, = a„ $± = ±|, sin£ = ±l. (45) 
We insert the ansatz (44) in Eq. (26) yielding 

^ s = 2irj s sech[2?7 s (v - Va )f] e -i[-p.(?-v.)t+*.+(**-°.)*]. ( 46 ) 

The fixed points of the phase portrait correspond to the time-independent i.e. v = v s . 
Using Eqs. (45), we obtain two fixed points 

Re^ = ±2r] s , Im^f = 0. (47) 

The ± signs refer to the two cases in Eqs. (45). 

Combining Eqs. (10) and (12) with Eq. (45), a s can be eliminated and we are left with 
two equations 

(27T 2 

v s = 2p s ± — — sechA s tanh A s , (48) 
-(K-p s )v s =p 2 s -4r ] 2 s -5± 

CITT 

— A s sechA s tanh A s , (49) 
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where A s = n(K— p s )/(4r) s ). For either sign, the system (48)-(49) has a single root v s = 2p s , 
p s = K and r/ s = — 5. Therefore, there are two fixed points on the real axis, located 

at Re^ = ±VK 2 - 5. 

For the set of parameters of Fig. 2 there are a stable and an unstable fixed point close 
to +1 and —1, respectively. The stability of the fixed points is determined by solving the 
collective coordinates equations numerically with the initial conditions very close to the 
above values, e.g., rj = r] s + 1CT 8 , q = 0, p = p s , $ = ®f = ±7r/2. In the unstable 
case ($o = — tt/2), the numerical solution exhibits oscillations of the amplitude and phase, 
whereas the velocity of the soliton remains constant. This solution is represented by the 
separatrix in Fig. 2. A trivial stable fixed point is located at the origin; its stability is 
established by numerical solutions of the collective coordinates equations with i] s close to 
zero. 
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